/*
 * File:   evolutionSP.cpp
 * Author: gian
 *
 * Created on December 11, 2012, 12:37 AM
 */

#include <math.h>
#include <iostream>
#include "evolutionSP.hpp"

evolutionSP::evolutionSP(double dh,
                         double dt) {
  mdh = dh;
  mdt = dt;
}

void evolutionSP::set_dt(double dt) {
  mdt = dt;
}

evolutionSP::evolutionSP(double dh,
                         double dt,
                         int    lsizex,
                         int    lsizey,
                         int    lsizez,
                         int    nx,
                         int    ny,
                         int    nz) {
  mdh              = dh;
  mdt              = dt;
  mLatticeSizeX    = lsizex;
  mLatticeSizeY    = lsizey;
  mLatticeSizeZ    = lsizez;
  mPartitionX      = nx;
  mPartitionY      = ny;
  mPartitionZ      = nz;
  mThreadNumber    = nx * ny * nz;
  mThreadArguments = new evolution_step_args[mThreadNumber];

  int deltax = mLatticeSizeX / nx;
  int deltay = mLatticeSizeY / ny;
  int deltaz = mLatticeSizeZ / nz;
  int i      = 0;

  for (int ix = 0; ix < nx; ix++) {
    for (int iy = 0; iy < ny; iy++) {
      for (int iz = 0; iz < nz; iz++) {
        mThreadArguments[i].evol = this;
        mThreadArguments[i].xmin = ix * deltax;
        mThreadArguments[i].xmax = (ix + 1) * deltax - 1;
        mThreadArguments[i].ymin = iy * deltay;
        mThreadArguments[i].ymax = (iy + 1) * deltay - 1;
        mThreadArguments[i].zmin = iz * deltaz;
        mThreadArguments[i].zmax = (iz + 1) * deltaz - 1;

        i++;
      }
    }
  }
}

evolutionSP::evolutionSP(const evolutionSP & orig) {}

evolutionSP::~evolutionSP() {
  delete[] mThreadArguments;
}

void evolutionSP::operator ()(latticeSP *          src,
                              latticeSP *          dst,
                              FriedmannEvolution * frw) {
  mLattice1 = src;
  mLattice2 = dst;
  
  mA0 = frw -> get_scale_parameter();
  frw->step(0.25*mdt);
  mA1 = frw -> get_scale_parameter();
  frw->step(0.5*mdt);    
  mA2 = frw -> get_scale_parameter();
  frw->step(0.25*mdt);
  
  for (int i = 0; i < mThreadNumber; i++) {
    pthread_create(mThreads + i, NULL, evolutionSP::evolution_step, mThreadArguments + i);
  }

  for (int i = 0; i < mThreadNumber; i++) {
    pthread_join(mThreads[i], NULL);
  }
}

void * evolutionSP::evolution_step(void * args) {
  evolution_step_args * p            = (evolution_step_args *) args;
  int                   lattice_size = p -> evol -> mLattice1 -> get_lattice_size();
  double                forces[p -> evol -> mLattice1 -> get_N()];
  double                force_tot, force_dominant;
  int                   k_dominant;
  int                   i, k;
  double                h         = p -> evol -> mdh;
  double                dt        = p -> evol -> mdt;
  double                a0        = p -> evol -> mA0;
  double                a1        = p -> evol -> mA1;
  double                a2        = p -> evol -> mA2;
  double                mul_pi    = dt*a0*h;
  double                mul_theta = 0.5*dt/(h*h*h)*(1.0/(a1*a1*a1)+1.0/(a2*a2*a2));

  
  for (int z = p -> zmin; z <= p -> zmax; z++) {
    for (int y = p -> ymin; y <= p -> ymax; y++) {
      for (int x = p -> xmin; x <= p -> xmax; x++) {
                    
        for (k = 0; k < p -> evol -> mLattice1 -> get_N(); k++) {
          forces[k] = 0.0;
        }

        force_tot = 0.0;

        nodeSP * psrc = p -> evol -> mLattice1 -> operator ()(x, y, z);

        // +x
        nodeSP * ptmp = p -> evol -> mLattice1 -> operator ()(x + 1, y, z);

        if (psrc -> N() == ptmp -> N()) {
          force_tot += sin(ptmp -> Theta() - psrc -> Theta());
        } else {
          forces[ptmp -> N()] -= sin(ptmp -> Theta() + psrc -> Theta());
        }

        // -x
        ptmp = p -> evol -> mLattice1 -> operator ()(x - 1, y, z);

        if (psrc -> N() == ptmp -> N()) {
          force_tot += sin(ptmp -> Theta() - psrc -> Theta());
        } else {
          forces[ptmp -> N()] -= sin(ptmp -> Theta() + psrc -> Theta());
        }

        // +y
        ptmp = p -> evol -> mLattice1 -> operator ()(x, y + 1, z);

        if (psrc -> N() == ptmp -> N()) {
          force_tot += sin(ptmp -> Theta() - psrc -> Theta());
        } else {
          forces[ptmp -> N()] -= sin(ptmp -> Theta() + psrc -> Theta());
        }

        // -y
        ptmp = p -> evol -> mLattice1 -> operator ()(x, y - 1, z);

        if (psrc -> N() == ptmp -> N()) {
          force_tot += sin(ptmp -> Theta() - psrc -> Theta());
        } else {
          forces[ptmp -> N()] -= sin(ptmp -> Theta() + psrc -> Theta());
        }

        // +z
        ptmp = p -> evol -> mLattice1 -> operator ()(x, y, z + 1);

        if (psrc -> N() == ptmp -> N()) {
          force_tot += sin(ptmp -> Theta() - psrc -> Theta());
        } else {
          forces[ptmp -> N()] -= sin(ptmp -> Theta() + psrc -> Theta());
        }

        // -z
        ptmp = p -> evol -> mLattice1 -> operator ()(x, y, z - 1);

        if (psrc -> N() == ptmp -> N()) {
          force_tot += sin(ptmp -> Theta() - psrc -> Theta());
        } else {
          forces[ptmp -> N()] -= sin(ptmp -> Theta() + psrc -> Theta());
        }

        for (k = 0; k < p -> evol -> mLattice1 -> get_N(); k++) {
          force_tot += forces[k];
        }

        double new_pi    = psrc -> Omega() + mul_pi * force_tot;
        double new_theta = psrc -> Theta() + mul_theta * new_pi;

        if (new_theta > M_PI) {
          new_theta      = 2.0 * M_PI - new_theta;
          new_pi         = -new_pi;
          force_dominant = forces[0];
          k_dominant     = 0;

          for (k = 1; k < p -> evol -> mLattice1 -> get_N(); k++) {
            if (forces[k] > force_dominant) {
              force_dominant = forces[k];
              k_dominant     = k;
            }
          }
        } else if (new_theta < 0.0) {
          new_theta      = -new_theta;
          new_pi         = -new_pi;
          force_dominant = forces[0];
          k_dominant     = 0;

          for (k = 1; k < p -> evol -> mLattice1 -> get_N(); k++) {
            if (forces[k] < force_dominant) {
              force_dominant = forces[k];
              k_dominant     = k;
            }
          }
        } else {
          k_dominant = psrc -> N();
        }

        nodeSP * pdst = p -> evol -> mLattice2 -> operator ()(x, y, z);

        pdst -> Omega() = new_pi;
        pdst -> Theta() = new_theta;
        pdst -> N()     = k_dominant;
      }
    }
  }

  pthread_exit(NULL);
}

double evolutionSP::deltaT() {
  return mdt;
}

double evolutionSP::deltaH() {
  return mdh;
}


//~ Formatted by Jindent --- http://www.jindent.com
